Scaling laws prediction from a solvable model 
of turbulent thermal convection 
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A solvable turbulent model is used to predict both the structure of the boundary layer and the 
scaling laws in thermal convection. The transport of heat depends on the interplay between the 
thermal, viscous and integral scales of turbulence, and thus, on both the Prandtl number and the 
Reynolds numbers. Depending on their values, a wide variety of possible regimes is found, including 
the classical 2/7 and 1/3 law, and a new 4/13 = 0.308 law for the Nusselt power law variation with 
the Rayleigh number. 

PACS numbers: 47.27 -i Turbulent flows, convection and heat transfer - 47.27.Eq Turbulence simulation and 
modeling - 47.27.Te Convection and heat transfer 
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Until the advent of high precision numerical and exper- 
imental data, the heat transport in turbulent convection 
was thought to be governed by the classical 1/3 power 
law linking the Nusselt number and the Rayleigh number. 
Deviations from this prediction have now been measured 
at large Rayleigh numbers, revealing a whole zoology of 
scaling exponents. They range from 0.2 to 0.25 in low 
Prandtl number, Mercury experiments [l|, to 0.29 ± 0.1 
in Helium experiments [2-4] , while electro-chemical con- 
vection gives the classical 1 /3 exponent Q ■ A number of 
theoretical models based on dimensional arguments have 
been developed to explain these new regimes [l|,[2|, [6fl, in- 
cluding the possibility that no real scaling prevails pj . In 
this letter, predictions obtained using a solvable model of 
turbulent convection are presented. This model couples 
large scale mean sheared velocity and temperature fields 
U = (£/(z),0,0) and 0(z), with small scale random ve- 
locity and temperature fields. This kind of large scale 
geometry is generally accepted as representative of the 
Boussinesq convection in the boundary layers, due to the 
shearing effect of the large convective cells. The model 
is closed by deriving an equation for the random compo- 
nent from the Boussinesq equation using two simplifying 
assumptions: i) the non-linear interactions of the small 
scale scales between themselves is modeled via a turbu- 
lent viscosity; ii) the small scale generation via the break- 
ing of large scale structures is modeled by a random small 
scale forcing with prescribed statistics. This results in a 
linear stochastic equation for the random small scales, 
which can be analytically solved in the shear flow geom- 
etry by a decomposition of the small waves into localized 
wave-packets. This model was used to obtain analytical 
predictions of mean velocity and/or temperature profiles 
in neutral boundary layer or channel flows [|[ or in the 
Planetary Surface Layer (PSL) @. Here, we adapt and 
generalize these results to determine both the structure 
of the boundary layers, and various scaling laws relating 
the Nusselt number Nu and various length scales to the 
Rayleigh number Ra. Some of our findings are then com- 
pared with results from high resolution direct numerical 
simulations (DNS) at Ra between 10 4 and 10 8 , and at 



Pr between 0.02 and 7, which are described in [10h12| . 

We consider Boussinesq equations, nondimensionalized 
by the thermal diffusivity and the cell height: 

d t Ui + dj {utUj) = -d L p + RaPr6S l3 + PrAm + ff' , 
d t e + d 3 {uj6) = A0 + /W. (1) 

Here, Pr is the Prandtl number, and and are 
small scale random forces which are introduced to model 
the seeding of small scales by turbulent plumes detach- 
ing from the wall and penetrating the outer region. For 
simplicity, this forcing is taken as spatially homogeneous 
and delta correlated in time (l3| . These assumptions 
only influence the magnitude of the Reynolds stresses, 
not their shape Q . Averaging ([T]) over the statistics 
of the forcing and assuming a shear flow geometry, one 
obtains the standard equations for the x-component of 
U =< u >= (U{z),0, 0) and for &(z) =< 9 >: 



d t U + d z < u'w' > = -d x P 
d t Q + d z < w'9' > = d 2 z O. 



PrdiU, 



(2) 



Here, the primes denote fluctuating (small scale) quan- 
tities and <> the averaging. In the stationary case, we 
get from ^ that d x P is a constant, independent of z. In 



w'9' > are neg- 



ligible, we thus obtain the well known parabolic profile 
for the velocity and the linear profile for the tempera- 
ture. In the turbulent case, the profiles are linear within 
the thermal or the viscous layer, while outside this layer, 
they are given by the condition 



d z < u'w' >= -d x P, d z < w'6' >= 0. 



(3) 



To close the system, we need < u'w' > and < w'6' >. 
For this, we now derive an equation for the fluctuating 
quantities, by taking into account the scale separation 
between the mean field and the random field l/L = r\ <C 
1. For this, we decompose the velocity field into localized 
wave-packets via a Gabor transform (GT): 



u(x,k,t) = y.g(?7*|x-x'|)e lk - (x - x ' ) u(x',t)dx', 



(4) 
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where g is a function which decreases rapidly at infinity 
and 1 > ?;* > r). Because of this, the GT of any large 
scale field is exponentially small in ry*, and the GT of 
Ujdju'i can be developed into Ujdju'i — dj^Umk^dkjU'i 
(see § for details). Using these properties, we obtain 
after GT of (fTJ) the small scale equation: 

(5) 



D t iii — —ikip — wd z USn + RaPr05i3 
D t 6 = -wd z Q -k 2 9 + p\ 



We have dropped primes on fluctuating quantities for 
convenient notations and introduced the total derivative 
Dt = dt + Ud x — d z (Uk x )dk z ■ In (0, we have further- 
more model the GT of the non-linear terms describing 
local interactions between fluctuations via the introduc- 
tion of turbulent viscosity and diffusivity, or equivalently, 
via a turbulent Prandtl number Pr*. The linear part of 
§5§ is exact and describes the non-local interactions be- 
tween the mean and the fluctuating part. We have thus 
obtained a linear, stochastic equation for the fluctuating 
quantities to close our system. 

To solve the closed system ^ and (O, we introduce 
two further simplification: i) we restrict ourselves to the 
2D case. This is justified because numerical work on 2D 
thermal convection showed that 2D geometry is sufficient 
to capture the ph ysic al mechanism responsible for the Nu 
scaling with Ra [L2\ . Analytical work on neutral shear 
flow has also proved that 2D geometry is sufficient to cap- 
ture the correct shape of the equilibrium profile [1] . Note 
that the vortex stretching, which is theoretically absent 
in 2D geometry, has been implicitly accounted for via the 
turbulent viscosity. In real 2D flows, this turbulent vis- 
cosity can be ignored [8j. ii) we take Pr* = I. This is 
because our model introduces a turbulent Prandtl num- 
ber whose exact value is unknown, which implies that the 
correct Pr scaling cannot be captured within our model. 
General Pr effects will be introduced in another manner, 
via links between length scales. With these approxima- 
tions, it was shown in 9] that to leading order in the 
Reynolds number Re, the solutions of ([5]) satisfy: 
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-(2+4a+)/3 



d z U 



2d z U(a + -a_) 2 V k 2 



' o x RaPr , f RaPr\ x 
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2d z U(a + -a_) 



d z U 
kl 



-(2+4a+)/3 



is the Richardson number. In developed convective tur- 
bulence, this number is large and negative R4 -C — 1. So 
only the leading order contribution in Ri in the expres- 
sion will be further considered. Combining ([3]) with 
the large Ri expansion of ©J we find that the turbulent 
profiles are given by the conditions: 



cx z, 



d z U 
u*V=R 
RaPr 



cx cte = Nu, 



(8) 



where we have introduced w 2 = (d z U /k 2 ) 2 ^~ Ri ^ . Look- 
ing for simple solutions where d z U and d z Q are power 
laws of z, the only solution is: 



d z U 



o z e 



1 



(9) 



resulting in a constant (with z) Richardson number. This 
solution corresponds to the standard logarithmic veloc- 
ity profile, and a temperature profile decreasin g lik e 1/z. 
Such a profile had been predicted by Malkus [14| using 
a maximum principle. It was found to be in good agree- 
ment with experiments of large Ra convection [l5j and 
DNS at Ra = I0 8 and Pr = 0.7 [U El- For the velocity 
profile, the prediction ([9]) is difficult to check numerically, 
because the Reynolds number is too low for the boundary 
layer to be fully develop ed flal . Even in the large Raleigh 
number experiments of [171 ] , the extent of the turbulent 
boundary layer is too small (a tenth of the cell) to check 
this prediction. However, in the PSL, the prediction © 
appears compatible with the measurements Q. 

It is also possible to find the leading order behavior in 
Re and Ri of vertical velocities and temperature fluctu- 
ations, using the results of @. They are: 



< w' 2 >« Ri 



d z U' 



< 



ld z U 



(RaPr) 2 ' 



(10) 



d z U 



8a+a 2 _—^——iJ,i + 4a+a_/X2 — 2a_ 
RaPr 



RaPr 

d r U 



This give a r.m.s. vertical velocity and temperature vary- 
ing like z 1 / 2 and z -1 / 2 , like in the free-fall regime. For the 
temperature fluctuation, the predicted decrease is more 
rapid than the classical z" 1 / 3 prediction. This feature 
has been already observed in high Ra convection [l7| ■ In 
the PSL, the r.m.s. velocity and temperature are com- 
patible with (Unj) @. 

From ((5]) and ©, we may also derive interesting exact 
relations. By matching the turbulent profiles with the 
viscous or diffusive solutions Prd z U = u 2 and d z Q = 
(JMj), where u T is the friction velocity and Nu the Nusselt 
mimber, we get: 



Here, k* is a characteristic horizontal wavenumber, Ai 
and Hi (i = 1, ..3) are some constant which depend 
on the forcing correlation functions, and a± = — (1 ± 
VI -4P,)/2, where 

R t = RaPr 



d z U 



2 Ay 

l T n > 

zPr 



d z e = Nu 



At 

z 



(11) 



o z e 
(d z uy 



(7) 



Here, we have introduced the thermal length scale At = 
1/Nu and the viscous length scale Ay = Pr/u T . To find 
the prefactor in (|5]), we use the law of energy dissipation 
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in a boundary layer geometry which gives e ~< u'w' > 
d z U. We then take into account the exact relation e = 
RaNu and use (jHJ, to obtain u 2 = RaNu and: 



< uw > 



RaPr 



Ra 



<w'6'> d z Uy/=Kl d z U 



(12) 



The first equality comes from (JSj> , while the second comes 
from the link between e and < u'w >. Condition (|12[) 
shows that Ri is independent of Ra. This condition is 
actually a necessary condition for stability of the turbu- 
lent boundary layer Q . The definition ([7]) combined with 
([TTj) then implies: 



Ra = uiNuPr. 



(13) 



We may also use (fTTj) and the above estimate to define 
a characteristic vertical velocity w c — y < w' 2 > / z and 
temperature fluctuation A c 
as: 



V< 0' 2 > z, which scale 



RaNu 



A,. 



u T Nu 
Ra 



(14) 



These exact relations derived within our nonlocal model 
will be the basis of the scaling theory which we now de- 
tail. 

For this, we follow [11| and introduce a third length 
scale, noted A/, representing the location of the peak 
of the kinetic horizontal energy spectrum. Its expres- 
sion can again be found within the nonlocal framework 
by noting that this length scale coincides with the loca- 
tion of the maximum of the r.m.s. horizontal velocity. 
In the non-local model, this horizontal velocity is pas- 
sively advected by the large scale horizontal velocity and 
obeys Ud x v! = d 2 u'. The solution to this equation de- 
pends on the structure of the velocity boundary layer: if 
the Reynolds number is too small, the boundary layer is 
mainly laminar, U — u r zj\y and simple dimensional ar- 
gument show that u' is a function of z/(x\vu^ 1 ) 1 ^ 3 . For 
a fixed aspect ratio T, this defines a typical vertical scale 
of r.m.s. velocity variation A/ = (TAyu" 1 ) 1 / 3 . When the 
Reynolds number exceeds a critical value of the order of 
Re c ~ 10 5 the boundary layer turbulence is well de- 
veloped, and most of the transport of vl is provided by 
the turbulent logarithmic regime U ~ u T [lnz + B). This 
regime is hard to obtain at Prandtl numbers of the order 
one or larger and it is likely that only low Prandtl num- 
ber system (like Mercury experiments) have the ability 
to reach this critical Reynolds number |l| . In this regime, 
u' is a function of z/jxu^ 1 ) 1 / 2 , thereby defining a typical 
length scale A/ ~ (Yu^: 1 ) 1 / 2 [19]. The two regimes can 
be lumped into the single formula A/ ~ Ay (Tu T Pr^ 2 )^ , 
where ip = 1/3 for Re < Re c and -0 = 1/2 for Re > Re c . 
Using these results, the logic of our argumentation is now 
to link u T to Nu by relating Ay and/or Aj to At, and 
then use the exact relation (fT3"j) to obtain Nu versus Ra. 
For this, we set At = \ v \)~ a , where < a < 1 is a 
parameter which enables to single out three remarkable 



regimes. In the first one, a = 1, At = Ay: the viscous 
and thermal layer coincide. This situation might be typi- 
cal of large aspect ratio cells, or in numerical simulations 
with stress-free boundary conditions 20]. In our model, 
it corresponds to Nu ~ Ra 1 / 3 Pr" 1 , the classical case. 
In the second regime, < a < 1 and the length-scale 
ordering changes to Ay < At < A/ . This situation seems 
to be typical for convection at Pr ~ 1 [11]. For illus- 
tration purpose, it is interesting to single out two special 
value of a which are relevant in turbulence: the first one 
is a = 2/3. This corresponds to At being the Taylor 
micro-scale, given by the square root of the mean energy 
to the mean enstrophy. The second value is a = 2/5. The 
corresponding scale varies like A/(A/i?e)- 3 / 10 , where Re 
is the Reynolds number based on the cell size and on the 
velocity at the integral scale Uj. It was shown [2l[ to 
correspond to the maximum of the Kolmogorov function 
and represents the location of the middle of the inertial 
range. At last, the third regime is obtained with a = 0, 
making At = A/ . This situation is typical of low Prandtl 
number convection. In our estimate of a, we have fol- 
lowed standard turbulence phenomenology, and assumed 
that A / = Ay (A/i?e) 3 / 4 . Note that for the heat trans- 
port, the exact relation in the two regimes ("0 = 1/3 or 
1/2) is: 



l-iHi-o 

Nu ~ Ra 3 -^ 1 -" 



2i/.[l-a] 5*I>\1 
- 3 _^|l_ Q j p r 3 



(15) 



The Prandtl dependence obtained in our model is 
stronger than in standard models of turbulent convec- 
tion. It is not necessarily inconsistent with available data. 
Such a dependence would for example account for most 
of the difference between the Pr = 4 experiment of [26[ 
and the Pr = 0.8 experiment of [22l |. However, we do not 
expect our model to fully capture the general Prandtl de- 
pendence, because of our approximation on the turbulent 
Prandtl number. In the sequel, we thus concentrate on 
the Rayleigh dependence of the physical quantities. The 
various scaling exponents predicted by the combination 
of (fT5|) , (fT4")l and the length scale relation is summarized 
in Table I. For purpose of comparison with the DNS, we 
have also included the predictions for Re < Re c at three 
special values of a, relevant to the low and order unity 
Prandtl number. 

In the low Prandtl regime, our model predicts a tran- 
sition between a Nu — Ra 1 / 4 regime up to Nu — Ra 1 / 5 
regime at larger Re or Ra. This is in agreement with the 
experimental findings of [l|. In the regime of Pr ~ 1, 
the prediction depends on the value of a. At Re < Re c , 
we find /3 = 4/13 (a = 2/3) and /3 = 2/7 (a = 2/5). 
The two values are close to the 0.29 value usually ob- 
served in experiments or in simulations. The 2/7 value 
is here exactly recovered when the thermal length scale 
matches the middle of the inertial range, which may help 
give a novel interpretation of the 2/7 law. The value 
j3 = 4/13 = 0.3077 is in remarkable agreement with the 
experimental value j3 = 0.309 ±0.0043 obtained in Q us- 
ing an Helium experiment with 0.5 aspect ratio and span- 
ning eleven decades of Rayleigh number. The case where 
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Name 



Nu 
Lw c /i> 
A c 
u T 
\v 
At 
Ai 



General case 



3-V>(l— a) 

1 3-24>(l-a) 

2 3-i/>(l-a) 

1 1 

2 3-i/>(l-a) 

1 

3-V(l-a) 

1 

3-l£(l-a) 
l-Vi(l-a) 
3-i/i(l-a) 
V— i) 

3-^(l-a) 



Value for i?e < i?e c = 1/3) 


a = 


a = 2/3 


a = 2/5 


0.2500 


0.3077 


0.2857 


0.4375 


0.4808 


0.4643 


-0.1875 


-0.1731 


-0.1786 


0.3750 


0.3462 


0.3571 


-0.3750 


-0.3462 


0.3571 


-0.2500 


-0.3077 


-0.2857 


-0.2500 


-0.2308 


-0.2381 



TABLE I: Summary of exponents as a function of the 
Rayleigh number as predicted by the present model. Here, 
tp = 1/3 when Re < Re c and ip = 1/2 when Re > Re c . 
< a < 1 is a parameter describing how the thermal 
length-scale compares with the viscous Ay or integral scale 
Xi, thereby describing different Prandtl regimes, like a — 
for low Prandtl, and a = 2/3 or 2/5 for Pr ~ 1. w c and A c 
are characteristic r.m.s. vertical velocity and temperature. 



our model. For w c and A c scaling with Ra, very good 
agreement between the prediction and the experimental 
value 0.43 and —0.2 [l| is obtained at low Prandtl num- 
ber. For Pr ~ 1, our estimates of w c L/v coincides for 
a = 2/3 with the value 0.485±0.005 measured in Helium 
0, while the exponent for A c is slightly larger than the 
experimental value —0.147. 

Another test of our predictions can also be made by 
comparison with DNS, where some length scales have 
been measured. The results are shown in Table II. They 
are in rather good agreement with the theoretical result, 
except for the case Pr = 7 of [1JJ where both the non- 
dimensional vertical velocity and the integral scale expo- 
nent deviate substantially from the prediction. It would 
certainly be interesting to investigate further the mean- 
ing of this discrepancy. For instance, it might be due to 
a Rayleigh dependence of the turbulent Prandtl number. 
Finally, we note that the aspect ratio is also likely to 
modify the recirculation pattern within a given cell and 



Re > Re c gives slightly different values. They are respec- 
tively (3 = 5/17 = 0.294 (a = 2/3) and /3 = 7/27 = 0.259 
(a = 2/5). The first value is so close to the value ob- 
tained before the transition, that it would make any tran- 
sition difficult to detect experimentally. This may explain 
why no transition was detected in [J], despite the large 
Rayleigh values attained. The second value is sufficiently 
lower than the value before transition, so that it should 
be detectable in a careful experiment. Finally, in the 
third regime, we find /? = 1/3 at all Reynolds number, 
in agreement with the Pr = 1 DNS of [l^j. Our model 
seems to rule out the (3 = 1/2 which has been predicted 
by Kraichnan 24] at very large Rayleigh number, and 
which might have been experimentally detected by Cha- 
vanne et al 25] . If this regime is confirmed, it might then 
correspond to a situation where the boundary layer grows 
unstable, thereby unvalidating a starting assumption of 



Name 


Verzicco and Camussi 


Werne (2D) 


Kerr and Herring 




Pr = 0.02 


Pr = 0.7 


Pr = 7 


Pr = 0.07 


Pr = 7 


Lw c /v 






0.54 


0.47 


0.56 


U T 






0.39 






At 


-0.25 


-0.29 




-0.26 


-0.29 


\i 


-0.18 


-0.23 




-0.26 


-0.11 



TABLE II: Summary of exponents as a function of the 
Rayleigh number as measured in various DNS at aspect ratio 
1 for the three first cases, and 4 for the last two. 



so different sca ling may appear at different aspect ratio, 
as observed in [261 ] . 

I thank Robert Kerr and Joseph Werne for references 
and interesting discussions. I have been partially funded 
by a NATO fellowship. 



[1] S. Cioni, S. Ciliberto and J. Sommeria J. Fluid Mech. 

335 111 (1997) 
[2] B. Castaing, G. Gunaratne, F. Heslot, L. Kadanoff, A. 

Libchaber,S. Thomae,X-Z. Wu, S. Zaleski and G. Zanetti 

J. Fluid Mech. 204 1 (1989). 
[3] T. Takeshita, T. Segawa, J.A. Glazier and M. Sano Phys. 

Rev. Letter 76 1465 (1996). 
[4] J.J. Niemela, L. Skrbek, K.R. Sreenivasan and R.J. Don- 
nelly preprint (1999). 
[5] R.J. Goldstein, HD. Chiang and D.L. See, J. Fluid Mech. 

213, 111 (1990) 
[6] B. I. Shraiman and ED. Siggia Phys. Rev. A. 42, 3650 

(1990). 

[7] S. Grossmann and D. Lohse J. FJuid Mech. 407 27 
(2000). 

[8] S. Nazarenko Phys. Lett. A 264 444 (2000); S. 
Nazarenko, N. Kevlahan and B. Dubrulle Physica D 139 
158 (2000); B. Dubrulle, J-P. Laval, S. Nazarenko and N. 



Kevlahan preprint. 
[9] B. Dubrulle, J-P. Laval and P. Sullivan preprint available 
at http://webast.ast.obs-mip.fr/people/bdubru (1999). 

[10] R. Verzicco and R. Camussi J. Fluid Mech. 383 55 (1999) 

[11] R.M. Kerr, J. Fluid Mech. 310 139 (1996); R. Kerr and 
J. Herring, submitted to J. Fluid Mech. 1999. 

[12] J. Werne, Phys. Rev. E 48 1020-1035 (1993). 

[13] Similar results would be obtained using a shortly corre- 
lated forcing. 

[14] W.V.R. Malkus Proc. Roy. Soc. London., A 225, 185- 

195. (1954); W.V.R. Malkus Proc. Roy. Soc. London., A 

225, 195-212. (1954) 
[15] A. A. Townsend, J. FJuid Mech. 5 209-241 (1959). 
[16] R. Kerr, private communication (1999). 
[17] A. Tilgner, A. Belmonte and A. Libchaber Phys. Rev 

E 47, R2253 (1993); A. Belmonte, A. Tilgner and A. 

Libchaber, Phys. Rev E 50 269 (1994). 
[18] H. Schlichting Boundary layer theory. Mc Graw Hill. 



5 



(1968) 

[19] A similar line of argumentation was done by pi for At, 

using the transport of O instead of it' . 
[20] J. Werne, private communication; K. Julien, S. Legg, 

J. Mc Williams and J. Werne, J. Fluid Mech. 322 243 

(1996). 

[21] F. Moisy, P. Tabeling and H. Willaime, Phys. Rev. Letter 

82 3994 (1999). 
[22] S. Ashkenazi and V. Steinberg Phys. Rev. Letter 83, 3641 

(1999). 



[23] X. Xu, K.M.S. Bajaj and G. Ahlers (2000), preprint. Such 
a dependence would for example account for most of the 
difference between the Pr = 4 experiment of [2(| and the 
Pr = 0.8 experiment of [H]. 
[24] R. Kraichnan Phys. Fluids 5 1374 (1962) 
[25] X. Chavanne, F. Chilla, B. Castaing, B. Hebral, B. 
Chabaud and J. Chaussy, Phys. Rev. Letters 79, 3648 
(1997). 

[26] X. Xu, K.M.S. Bajaj and G. Ahlers (2000), preprint. 



